#include "StdAfx.h"
#include "Kmeans.h"
#include <complex>
#define natureNum 2.7182818284590452353602874713527

Kmeans::Kmeans(int inputD, double** input, int  NumofCluster, int inputS)
{
	inputSize = inputS;
	inputData = input;
	inputDim = inputD;
	numCluster = NumofCluster;

	Kmeans::norInput();

	numIteration = 0;

	newSolution = (int*)malloc(sizeof(int)*inputSize);
	currentSolution = (int*)malloc(sizeof(int)*inputSize);
	finalSolution = (int*)malloc(sizeof(int)*inputSize);

	currentCentroid = (double**)malloc(sizeof(double*)*numCluster);
	newCentroid = (double**)malloc(sizeof(double*)*numCluster);
	finalCentroid = (double**)malloc(sizeof(double*)*numCluster);
	currentSum = (double**)malloc(sizeof(double*)*numCluster);
	newSum = (double**)malloc(sizeof(double*)*numCluster);
	currentElement = (int*)calloc(numCluster, sizeof(int));
	newElement = (int*)calloc(numCluster, sizeof(int));
	srand((unsigned)time(NULL));
	for(int i=0; i<numCluster; i++)
	{
		currentCentroid[i] = (double*)malloc(sizeof(double)*inputDim);
		newCentroid[i] = (double*)malloc(sizeof(double)*inputDim);
		finalCentroid[i] = (double*)malloc(sizeof(double)*inputDim);
		currentSum[i] = (double*)calloc(inputDim, sizeof(double));
		newSum[i] = (double*)calloc(inputDim, sizeof(double));
	}

	Kmeans::findSolution();
	finalCost = currentCost;
	for(int i=0; i<inputSize; i++)
	{
		finalSolution[i] = currentSolution[i];
	}
	for(int i=0; i<numCluster; i++)
	{
		for(int j=0; j<inputDim; j++)
		{
			finalCentroid[i][j] = currentCentroid[i][j];
		}
	}
	for(int k=0; k<32; k++)
	{
		Kmeans::findSolution();
		if(finalCost<currentCost)
		{
			finalCost = currentCost;
			for(int i=0; i<inputSize; i++)
			{
				finalSolution[i] = currentSolution[i];
			}
			for(int i=0; i<numCluster; i++)
			{
				for(int j=0; j<inputDim; j++)
				{
					finalCentroid[i][j] = currentCentroid[i][j];
				}
			}
		}
	}
/*	CString str;
	for(int i=0; i<inputSize; i++)
	{
		str += douTostr(inputData[i][0]);
		str += "\t";
		str += douTostr(finalSolution[i]);
		str += "\n";
	}
	str += "\n";
	str += douTostr(finalCost);
	str += "\t";
	str += douTostr(numIteration);
	::AfxMessageBox(str, MB_OK, 0);*/
}
//initial a solution by random assign
void Kmeans::initialSolution(void)
{
	for(int i=0; i<inputSize; i++)
	{
		currentSolution[i] = rand()%numCluster;
		newSolution[i] = currentSolution[i];
	}
	//get the centroid
	for(int i=0; i<numCluster; i++)
	{
		currentElement[i] = 0;
		newElement[i] = 0;
		for(int j=0; j<inputDim; j++)
		{
			currentSum[i][j] = 0;
			newSum[i][j] = 0;
		}
	}

	for(int i=0; i<inputSize; i++)
	{
		currentElement[currentSolution[i]] ++;
		newElement[currentSolution[i]] ++;
		for(int j=0; j<inputDim; j++)
		{
			currentSum[currentSolution[i]][j] += inputData[i][j];
			newSum[currentSolution[i]][j] += inputData[i][j];
		}
	}

	for(int i=0; i<numCluster; i++)
	{
		for(int j=0; j<inputDim; j++)
		{
			currentCentroid[i][j] = currentSum[i][j]/((currentElement[i]==0)?1:currentElement[i]);
			newCentroid[i][j] = newSum[i][j]/((newElement[i]==0)?1:newElement[i]);
		}
	}
	return;
}
//search for the optimal solution
void Kmeans::findSolution(void)
{
	Kmeans::initialSolution(); //initial clustering
	currentCost = Kmeans::getInitialCost();
	//recursive update
	bool flag = 1;
	while(flag)
	{
		numIteration ++;
		Kmeans::updateSolution();
		newCost = Kmeans::getNewCost();
		if(newCost<currentCost)
		{
			currentCost = newCost;
			for(int i=0; i<inputSize; i++)
			{
				currentSolution[i] = newSolution[i];
			}
			for(int i=0; i<numCluster; i++)
			{
				currentElement[i] = newElement[i];
				for(int j=0; j<inputDim; j++)
				{
					currentSum[i][j] = newSum[i][j];
					currentCentroid[i][j] = newCentroid[i][j];
				}
			}
		}
		else
			flag = 0;
	}
	//initial annealing parameters
/*	currentT = 1;
	coolingRate = 0.996;
	minimumT = 0.05;
	//start annealing
	while(currentT>minimumT)
	{
		for(int kk=0; kk<5; kk++)
		{
		numIteration ++;
		Kmeans::sAnnealing(); //get a new state
		newCost = Kmeans::getNewCost(); //get the new cost
		double changedCost = newCost - currentCost; //compute the delta cost
		if(numIteration%1 == 0)
		{
			CString str;
			for(int i=0; i<inputSize; i++)
			{
				str += douTostr(inputData[i][0]);
				str += "\t";
				str += douTostr(currentSolution[i]);
				str += "\t";
				str += douTostr(newSolution[i]);
				str += "\n";
			}
			str += douTostr(currentT);
			str += "\t";
			str += douTostr(currentCost);
			str += "\t";
			str += douTostr(newCost);
			str += "\n";
			str += douTostr(numIteration);
			::AfxMessageBox(str, MB_OK, 0);
		}		if(changedCost<=0) 
		{
			Kmeans::takeChange();
		}
		else
		{
			double temp = pow(natureNum, -changedCost/currentT);
			double tempy = (double)((unsigned)rand())/(double)((unsigned)rand());
			tempy = (tempy<1)?temp:(1/tempy); 
			if(temp>tempy)
			{
				Kmeans::takeChange();
			}
			else
			{
				Kmeans::rejectChange();
			}
		}
		}
		//change the parameters
		currentT *= coolingRate;
	}*/
	return;
}
//take the change
void Kmeans::takeChange(void)
{
	currentSolution[changedInput] = newCluster;
	currentElement[oldCluster] = newElement[oldCluster];
	currentElement[newCluster] = newElement[newCluster];
	for(int j=0; j<inputDim; j++)
	{
		currentSum[oldCluster][j] = newSum[oldCluster][j];
		currentSum[newCluster][j] = newSum[oldCluster][j];
		currentCentroid[oldCluster][j] = newCentroid[oldCluster][j];
		currentCentroid[newCluster][j] = newCentroid[oldCluster][j];
	}
	currentCost = newCost;
}
//reject the change
void Kmeans::rejectChange(void)
{
	newSolution[changedInput] = oldCluster;
	newElement[oldCluster] = currentElement[oldCluster];
	newElement[newCluster] = currentElement[newCluster];
	for(int j=0; j<inputDim; j++)
	{
		newSum[oldCluster][j] = currentSum[oldCluster][j];
		newSum[oldCluster][j] = currentSum[newCluster][j];
		newCentroid[oldCluster][j] = currentCentroid[oldCluster][j];
		newCentroid[newCluster][j] = currentCentroid[newCluster][j];
	}
}
//update the solution
void Kmeans::updateSolution(void)
{
	for(int i=0; i<inputSize; i++)
	{
		newSolution[i] = Kmeans::findCluster(inputData[i]);
	}
	//get the centroid
	for(int i=0; i<numCluster; i++)
	{
		newElement[i] = 0;
		for(int j=0; j<inputDim; j++)
		{
			newSum[i][j] = 0;
		}
	}

	for(int i=0; i<inputSize; i++)
	{
		newElement[newSolution[i]] ++;
		for(int j=0; j<inputDim; j++)
		{
			newSum[newSolution[i]][j] += inputData[i][j];
		}
	}

	for(int i=0; i<numCluster; i++)
	{
		for(int j=0; j<inputDim; j++)
		{
			newCentroid[i][j] = newSum[i][j]/((newElement[i]==0)?1:newElement[i]);
		}
	}
	return;
}
//find the proper cluster
int Kmeans::findCluster(double* currentinput)
{
	int currentCluster;
	double minDist = 999999999;
	for(int i=0; i<numCluster; i++)
	{
		double temp = Kmeans::wDistance(currentinput, currentCentroid[i]);
		if(temp < minDist)
		{
			minDist = temp;
			currentCluster = i;
		}
	}
	return currentCluster;
}
//annleaing
void Kmeans::sAnnealing(void)
{
	changedInput = rand()%inputSize;
	newCluster = rand()%numCluster;
	oldCluster = newSolution[changedInput];
	newSolution[changedInput] = newCluster;
	newElement[oldCluster] --;
	newElement[newCluster] ++;
	for(int j=0; j<inputDim; j++)
	{
		newSum[oldCluster][j] -= inputData[changedInput][j];
		newSum[newCluster][j] += inputData[changedInput][j];
		newCentroid[oldCluster][j] = newSum[oldCluster][j] / ((newElement[oldCluster]==0)?1:newElement[oldCluster]);
		newCentroid[newCluster][j] = newSum[newCluster][j] / newElement[newCluster];
	}
}
//compute the distance
double Kmeans::wDistance(double* x, double* y)
{
	double disttemp = 0;
	for(int i=0; i<inputDim; i++)
	{
		disttemp += pow(x[i]-y[i], 2);
	}
	return sqrt(disttemp);
}
//get the initial cost
double Kmeans::getInitialCost(void)
{
	double sum = 0;
	for(int i=0; i<inputSize; i++)
	{
		sum += Kmeans::wDistance(inputData[i], currentCentroid[currentSolution[i]]);
	}
	return sum;
}
//get the new cost
double Kmeans::getNewCost(void)
{
	double sum = 0;
	for(int i=0; i<inputSize; i++)
	{
		sum += Kmeans::wDistance(inputData[i], newCentroid[newSolution[i]]);
	}
	return sum;
}
//normalize input
void Kmeans::norInput(void)
{
	double* inputSum = (double*)calloc(inputDim, sizeof(double));
	double* inputAve = (double*)malloc(sizeof(double)*inputDim);
	double* inputStd = (double*)malloc(sizeof(double)*inputDim);
	//get the sum
	for(int i=0; i<inputSize; i++)
	{
		for(int j=0; j<inputDim; j++)
		{
			inputSum[j] += inputData[i][j];
		}
	}
	//get the average
	for(int j=0; j<inputDim; j++)
	{
		inputAve[j] = inputSum[j]/inputSize;
		inputSum[j] = 0;
	}
	//calculate the std
	for(int i=0; i<inputSize; i++)
	{
		for(int j=0; j<inputDim; j++)
		{
			inputSum[j] += pow(inputData[i][j]-inputAve[j], 2);
		}
	}
	for(int j=0; j<inputDim; j++)
	{
		inputStd[j] = sqrt(inputSum[j]/(inputSize-1));
	}
	for(int j=0; j<inputDim; j++)
	{
		if(inputStd[j]!=0)
		{
			for(int i=0; i<inputSize; i++)
			{
				inputData[i][j] = (inputData[i][j]-inputAve[j])/inputStd[j];
			}
		}
	}
	free(inputSum);
	free(inputAve);
	free(inputStd);
	return;
}
//get the results
int* Kmeans::getResluts(void)
{
	return finalSolution;
}
//get the centroid
double** Kmeans::getCentroid(void)
{
	return finalCentroid;
}
//get the sum distance
double Kmeans::getCost(void)
{
	return finalCost;
}
//get the number of iteration
int Kmeans::getIteration(void)
{
	return numIteration;
}
//get the test results
int Kmeans::getTest(double* testData)
{
	double mindist = 99999999;
	int testR = 0;
	for(int i=0; i<numCluster; i++)
	{
		double temp = Kmeans::wDistance(testData, finalCentroid[i]);
		if(temp<mindist)
		{
			mindist = temp;
			testR = i;
		}
	}
	return testR;
}
Kmeans::~Kmeans(void)
{
}

//convert double to string
CString Kmeans::douTostr( double dou )
{
	int decimal , sign, i;
	char *buffer = new char [10];
	char *buffer1 = new char [10] ;
	int precision = 10;
	int err;

	err = _ecvt_s( buffer, _CVTBUFSIZE, dou , precision , &decimal , & sign ) ;

	if ( decimal == 0 )
	{
		buffer1[0] = '0' ;
		buffer1[1] = '.' ;
		for ( int i = 2 ; i < 10 ; i ++ )
		{
			buffer1[i] = buffer[i-2] ;
		}
	}
	if ( decimal > 0 )
	{
		for ( i = 0 ; i < decimal ; i ++ )
		{
			buffer1[i] = buffer[i] ;
		}
		buffer1[i++] = '.' ;
		for ( ; i < 10 ; i ++ )
		{
			buffer1[i] = buffer[i-1] ;
		}
	}
	if ( decimal < 0 )
	{
		buffer1[0] = '0' ;
		buffer1[1] = '.' ;
		for ( i = 2 ; i < abs ( decimal ) + 2 ; i ++ )
		{
			buffer1[i] = '0' ;
		}
		int j = 0 ;
		for ( ; i < 10 ; i ++ )
		{
			buffer1[i] = buffer[j++] ;
		}
	}

	CString cstring ;

	if(dou<0)
		cstring = "-";
	for ( i = 0 ; i < 10 ; i ++ )
	{
		cstring = cstring + buffer1[i] ;
	}
	return cstring ;
}
